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O ' Abstract 

(N ' 

Signals comprised of a stream of short pulses appear in many applications including bio-imaging and 

radar. The recent finite rate of innovation framework, has paved the way to low rate sampling of such 

pulses by noticing that only a small number of parameters per unit time are needed to fully describe 

these signals. Unfortunately, for high rates of innovation, existing sampling schemes are numerically 

unstable. In this paper we propose a general sampling approach which leads to stable recovery even in the 

presence of many pulses. We begin by deriving a condition on the sampling kernel which allows perfect 

reconstruction of periodic streams from the minimal number of samples. We then design a compactly 

| supported class of filters, satisfying this condition. The periodic solution is extended to finite and infinite 

CNJ . streams, and is shown to be numerically stable even for a large number of pulses. High noise robustness 

00 

■ is also demonstrated when the delays are sufficiently separated. Finally, we process ultrasound imaging 

data using our techniques, and show that substantial rate reduction with respect to traditional ultrasound 

o ; 

sampling schemes can be achieved. 

Index Terms 

X 

Analog-to-digital conversion, annihilating filters, finite rate of innovation, compressed sensing, perfect 
reconstruction, ultrasound imaging, sub-Nyquist sampling. 



Copyright (c) 2010 IEEE. Personal use of this material is permitted. However, permission to use this material for any other 
purposes must be obtained from the IEEE by sending a request to pubs-permissions@ieee.org. 

Department of Electrical Engineering, Technion — Israel Institute of Technology, Haifa 32000, Israel. Phone: +972-4-8293256, 
fax: +972-4-8295757, E-mail: {ronentur@techunix,yonina@ee}.technion.ac.il, zvi.friedman@med.ge.com. Y. Eldar is currently 
a Visiting Professor at Stanford, CA. This work was supported in part by a Magneton grant from the Israel Ministry of Industry 
and Trade. 



2 



I. Introduction 

Sampling is the process of representing a continuous-time signal by discrete-time coefficients, while 
retaining the important signal features. The well-known Shannon-Nyquist theorem states that the minimal 
sampling rate required for perfect reconstruction of bandlimited signals is twice the maximal frequency. 
This result has since been generalized to minimal rate sampling schemes for signals lying in arbitrary 
subspaces ifTI. El. 

Recently, there has been growing interest in sampling of signals consisting of a stream of short pulses, 
where the pulse shape is known. Such signals have a finite number of degrees of freedom per unit time, 
also known as the Finite Rate of Innovation (FRI) property Q- This interest is motivated by applications 
such as digital processing of neuronal signals, bio-imaging, image processing and ultrawideband (UWB) 
communications, where such signals are present in abundance. Our work is motivated by the possible 
application of this model in ultrasound imaging, where echoes of the transmit pulse are reflected off 
scatterers within the tissue, and form a stream of pulses signal at the receiver. The time-delays and am- 
plitudes of the echoes indicate the position and strength of the various scatterers, respectively. Therefore, 
determining these parameters from low rate samples of the received signal is an important problem. 
Reducing the rate allows more efficient processing which can translate to power and size reduction of 
the ultrasound imaging system. 

Our goal is to design a minimal rate single-channel sampling and reconstruction scheme for pulse 
streams that is stable even in the presence of many pulses. Since the set of FRI signals does not form 
a subspace, classic subspace schemes cannot be directly used to design low-rate sampling schemes. 
Mathematically, such FRI signals conform with a broader model of signals lying in a union of subspaces 
Although the minimal sampling rate required for such settings has been derived, no generic 
sampling scheme exists for the general problem. Nonetheless, some special cases have been treated in 
previous work, including streams of pulses. 

A stream of pulses can be viewed as a parametric signal, uniquely defined by the time-delays of the 
pulses and their amplitudes. Efficient sampling of periodic impulse streams, having L impulses in each 
period, was proposed in Q, iflOl . The heart of the solution is to obtain a set of Fourier series coefficients, 
which then converts the problem of determining the time-delays and amplitudes to that of finding the 
frequencies and amplitudes of a sum of sinusoids. The latter is a standard problem in spectral analysis 
ifTTI which can be solved using conventional methods, such as the annihilating filter approach, as long 
as the number of samples is at least 2L. This result is intuitive since there are 2L degrees of freedom in 
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each period: L time-delays and L amplitudes. 

Periodic streams of pulses are mathematically convenient to analyze, however not very practical. In 
contrast, finite streams of pulses are prevalent in applications such as ultrasound imaging. The first 
treatment of finite Dirac streams appears in O, in which a Gaussian sampling kernel was proposed. The 
time-delays and amplitudes are then estimated from the Gaussian tails. This method and its improvement 
|[T2l are numerically unstable for high rates of innovation, since they rely on the Gaussian tails which take 
on small values. The work in |[T3l introduced a general family of polynomial and exponential reproducing 
kernels, which can be used to solve FRI problems. Specifically, B-spline and E-spline sampling kernels 
which satisfy the reproduction condition are proposed. This method treats streams of Diracs, differentiated 
Diracs, and short pulses with compact support. However, the proposed sampling filters result in poor 
reconstruction results for large L. To the best of our knowledge, a numerically stable sampling and 
reconstruction scheme for high order problems has not yet been reported. 

Infinite streams of pulses arise in applications such as UWB communications, where the communicated 
data changes frequently. Using spline filters |[T3l , and under certain limitations on the signal, the infinite 
stream can be divided into a sequence of separate finite problems. The individual finite cases may be 
treated using methods for the finite setting, at the expense of above critical sampling rate, and suffer 
from the same instability issues. In addition, the constraints that are cast on the signal become more and 
more stringent as the number of pulses per unit time grows. In a recent work lfl4l the authors propose 
a sampling and reconstruction scheme for L = 1, however, our interest here is in high values of L. 

Another related work Q proposes a semi-periodic model, where the pulse time-delays do not change 
from period to period, but the amplitudes vary. This is a hybrid case in which the number of degrees 
of freedom in the time-delays is finite, but there is an infinite number of degrees of freedom in the 
amplitudes. Therefore, the proposed recovery scheme generally requires an infinite number of samples. 
This differs from the periodic and finite cases we discuss in this paper which have a finite number of 
degrees of freedom and, consequently, require only a finite number of samples. 

In this paper we study sampling of signals consisting of a stream of pulses, covering the three different 
cases: periodic, finite and infinite streams of pulses. The criteria we consider for designing such systems 
are: a) Minimal sampling rate which allows perfect reconstruction, b) numerical stability (with sufficiently 
separated time delays), and c) minimal restrictions on the number of pulses per sampling period. 

We begin by treating periodic pulse streams. For this setting, we develop a general sampling scheme 
for arbitrary pulse shapes which allows to determine the times and amplitudes of the pulses, from a 
minimal number of samples. As we show, previous work J3J is a special case of our extended results. 
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In contrast to the infinite time-support of the filters in Q, we develop a compactly supported class of 
filters which satisfy our mathematical condition. This class of filters consists of a sum of sine functions 
in the frequency domain. We therefore refer to such functions as Sum of Sines (SoS). To the best of our 
knowledge, this is the first class of finite support filters that solve the periodic case. As we discuss in 
detail in Section [V] these filters are related to exponential reproducing kernels, introduced in |[T3l . 

The compact support of the SoS filters is the key to extending the periodic solution to the finite stream 
case. Generalizing the SoS class, we design a sampling and reconstruction scheme which perfectly 
reconstructs a finite stream of pulses from a minimal number of samples, as long as the pulse shape has 
compact support. Our reconstruction is numerically stable for both small values of L and large number 
of pulses, e.g., L = 100. In contrast, Gaussian sampling filters [31, lfT2l are unstable for L > 9, and 
we show in simulations that B-splines and E-splines |fT3l exhibit large estimation errors for L > 5. 
In addition, we demonstrate substantial improvement in noise robustness even for low values of L. Our 
advantage stems from the fact that we propose compactly supported filters on the one hand, while staying 
within the regime of Fourier coefficients reconstruction on the other hand. Extending our results to the 
infinite setting, we consider an infinite stream consisting of pulse bursts, where each burst contains a 
large number of pulses. The stability of our method allows to reconstruct even a large number of closely 
spaced pulses, which cannot be treated using existing solutions ifTBl . In addition, the constraints cast 
on the structure of the signal are independent of L (the number of pulses in each burst), in contrast to 
previous work, and therefore similar sampling schemes may be used for different values of L. Finally, 
we show that our sampling scheme requires lower sampling rate for L > 3. 

As an application, we demonstrate our sampling scheme on real ultrasound imaging data acquired by 
GE healthcare's ultrasound system. We obtain high accuracy estimation while reducing the number of 
samples by two orders of magnitude in comparison with current imaging techniques. 

The remainder of the paper is organized as follows. In Section [TT] we present the periodic signal 
model, and derive a general sampling scheme. The SoS class is then developed and demonstrated via 
simulations. The extension to the finite case is presented in Section [TTTl followed by simulations showing 
the advantages of our method in high order problems and noisy settings. In Section [IV] we treat infinite 
streams of pulses. Section [V] explores the relationship of our work to previous methods. Finally, in 
Section [V]] we demonstrate our algorithm on real ultrasound imaging data. 
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II. Periodic Stream of Pulses 

A. Problem Formulation 

Throughout the paper we denote matrices and vectors by bold font, with lowercase letters corresponding 
to vectors and uppercase letters to matrices. The nth element of a vector a is written as a n , and Ajj 
denotes the ijth element of a matrix A. Superscripts (•)*, (■) T and (-) H represent complex conjugation, 
transposition and conjugate transposition, respectively. The Moore-Penrose pseudo-inverse of a matrix A 
is written as A^. The continuous-time Fourier transform (CTFT) of a continuous-time signal x (t) € L2 
is defined by X (u) = JZo x (*) e~ jujt dt, and 

/oo 
x*(t)y(t)dt, (1) 
-00 

denotes the inner product between two L2 signals. 

Consider a r-periodic stream of pulses, defined as 

L 

x(t) = ^ ^Ja^(t - - mr), (2) 

m£Z 1=1 

where h(t) is a known pulse shape, r is the known period, and {ti, ai}f =l , ti G [0,t), a\ G C, I = 1 . . . L 
are the unknown delays and amplitudes. Our goal is to sample x(t) and reconstruct it, from a minimal 
number of samples. Since the signal has 2L degrees of freedom, we expect the minimal number of 
samples to be 2L. We are primarily interested in pulses which have small time-support. Direct uniform 
sampling of 2L samples of the signal will result in many zero samples, since the probability for the 
sample to hit a pulse is very low. Therefore, we must construct a more sophisticated sampling scheme. 

Define the periodic continuation of h(t) as f(t) = ^ m£Z li(t — mr). Using Poisson's summation 
formula |fl"5l , f(t) may be written as 

where H(cu) denotes the CTFT of the pulse h(t). Substituting ([3]) into Q we obtain 

L 

i=i 




fcez 
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where we denoted 

X [k] = -H (—\ aie~^ kt ^. (5) 
r ^ T ' 1=1 

The expansion in (O is the Fourier series representation of the r-periodic signal x(t) with Fourier 
coefficients given by ((5]). 

Following 0, we now show that once 2L or more Fourier coefficients of x(t) are known, we may use 
conventional tools from spectral analysis to determine the unknowns {ti,ai}f =1 . The method by which 
the Fourier coefficients are obtained will be presented in subsequent sections. 

Define a set K, of M consecutive indices such that H (^f^) / 0, Vfc € /C. We assume such a set 
exists, which is usually the case for short time-support pulses h(t). Denote by H the M X M diagonal 
matrix with kth entry (— ), and by V(t) the M x L matrix with kith element e -'j 2 ^ kt i/ T j where 
t = {t\, . . . , it,} is the vector of the unknown delays. In addition denote by a the length-L vector whose 
/th element is ai, and by x the length-M vector whose fcfh element is X[k]. We may then write ([5]) in 
matrix form as 

x = HV(t)a. (6) 

Since H is invertible by construction we define y = H~ 1 x, which satisfies 

y = V(t)a. (7) 

The matrix V is a Vandermonde matrix and therefore has full column rank ifTTIl . |[T6l as long as M > L 
and the time-delays are distinct, i.e., ^ tj for all i ^ j. 

Writing the expression for the kth element of the vector y in ^} explicitly: 

L 

y fc = J]a z e-^/ r . (8) 
l=i 

Evidently, given the vector x, ((TJ is a standard problem of finding the frequencies and amplitudes of a 
sum of L complex exponentials (see ifTTI for a review of this topic). This problem may be solved as 
long as \K\ = M > 2L. 

The annihilating filter approach used extensively by Vetterli et al. 0, |[T0l is one way of recovering 
the frequencies, and is thoroughly described in the literature 0, ifTOl . IfTTI . This method can solve 
the problem using the critical number of samples M = 2L, as opposed to other techniques such as 
MUSIC Ifm , ifTBl and ESPRIT |[T9l which require oversampling. Since we are interested in minimal -rate 
sampling, we use the annihilating filter throughout the paper. 
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B. Obtaining The Fourier Series Coefficients 

As we have seen, given the vector of M > 2L Fourier series coefficients x, we may use standard 
tools from spectral analysis to determine the set {tua{\f =l . In practice, however, the signal is sampled 
in the time domain, and therefore we do not have direct access to samples of x. Our goal is to design a 
single-channel sampling scheme which allows to determine x from time-domain samples. In contrast to 
previous work 0> iTTOll which focused on a low-pass sampling filter, in this section we derive a general 
condition on the sampling kernel allowing to obtain the vector x. For the sake of clarity we confine 
ourselves to uniform sampling, the results extend in a straightforward manner to nonuniform sampling 
as well. 



x(t) 



s*(-t) 



t = nT 



Fig. 1. Single channel sampling scheme. 



Consider sampling the signal x(t) uniformly with sampling kernel s*(—t) and sampling period T, as 
depicted in Fig. Q] The samples are given by 



cm 



Substituting © into ® we have 



c n 



x(t)s*(t - nT)dt = (s(t - nT),x(t)). 



/oo 
e j27Tkt/T s*{t-nT)dt 

/OO 
e j2jTkt/T s*(t)dt 
-oo 

^X[k]e j2nknT/r S*{2Trk/T), 



where S(cj) is the CTFT of s(t). Choosing any filter s(t) which satisfies 

uj = 2-nk/r, k (£ K 



S(oj) = < 



nonzero u = 2-Kk/r, k G K, 
arbitrary otherwise, 



we can rewrite (fTD~b as 



c [n] = J2x[k]e j27TknT/r S*(2TTk/T). 



(9) 



(10) 



(11) 



(12) 
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In contrast to (fTOb . the sum in (fT2l ) is finite. Note that (fTTT) implies that any real filter meeting this 
condition will satisfy k G K, =^ — A; G /C, and in addition S(2irk/T) = S*(—2nk/T), due to the conjugate 
symmetry of real filters. 

Defining the M x M diagonal matrix S whose Mi entry is S*(2irk/T) for all k G /C, and the length- iV 
vector c whose nth element is c[n], we may write (fT2l ) as 

c = V(-t s )Sx (13) 

where t s = {nT : n = . . . N — 1}, and V is defined as in © with a different parameter — t s and 
dimensions N x M. The matrix S is invertible by construction. Since V is Vandermonde, it is left 
invertible as long as N > M. Therefore, 

x = S- 1 Vt(-t s )c. (14) 

In the special case where N = M and T = r/N, the recovery in (fl4l) becomes: 

x = S^DFT-tc}, (15) 

i.e., the vector x is obtained by applying the Discrete Fourier Transform (DFT) on the sample vector, 
followed by a correction matrix related to the sampling filter. 

The idea behind this sampling scheme is that each sample is actually a linear combination of the 
elements of x. The sampling kernel s(t) is designed to pass the coefficients X[k], k G K, while suppressing 
all other coefficients X[k], k ^ /C. This is exactly what the condition in (fTTb means. This sampling scheme 
guarantees that each sample combination is linearly independent of the others. Therefore, the linear system 
of equations in (fT3l has full column rank which allows to solve for the vector x. 

We summarize this result in the following theorem. 

Theorem 1. Consider the r-periodic stream of pulses of order L: 

L 

x(t) = aih(t — ti — mr). 

m£Z 1=1 

Choose a set K, of consecutive indices for which H(2irk/T) ^ 0, \/k G M Then the samples 

c[n] = (s(t-nT),x(t)), n = . . . N - 1, 
uniquely determine the signal x(t) for any s(t) satisfying condition (II II ). as long as N > \K\ > 2L. 

In order to extend Theorem Q] to nonuniform sampling, we only need to substitute the nonuniform 
sampling times in the vector t s in (fT4l ). 
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Fig. 2. The filter g(t) with all coefficients bk = 1. 

Theorem Q] presents a general single channel sampling scheme. One special case of this framework 
is the one proposed by Vetterli et al. in in which s*(—t) = B smc(-Bt), where B = M/t and 
N > M > 2L. In this case s(t) is an ideal low-pass filter of bandwidth B with 



Clearly, (fJU) satisfies the general condition in CCD with K = {-[M/2\,. . . , [M/2\} and S (^) = 
-7=, V/c G /C. Note that since this filter is real valued it must satisfy fce/C=^-— k £ JC, i.e., the indices 

V 2"7T 

come in pairs except for /c = 0. Since k = is part of the set /C, in this case the cardinality M = \K\ 
must be odd valued so that iV > M > 2L + 1 samples, rather than the minimal rate N > 2L. 

The ideal low-pass filter is bandlimited, and therefore has infinite time-support, so that it cannot be 
extended to finite and infinite streams of pulses. In the next section we propose a class of non-bandlimited 
sampling kernels, which exploit the additional degrees of freedom in condition (fTTT) . and have compact 
support in the time domain. The compact support allows to extend this class to finite and infinite streams, 
as we show in Sections [III] and [IV] respectively. 

C. Compactly Supported Sampling Kernels 

Consider the following SoS class which consists of a sum of sines in the frequency domain: 




(16) 




(17) 
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Fig. 3. The filter g(t) with Hamming window coefficients. 



where bt ^ 0, k G K. The filter in (fTTT ) is real valued if and only if k G K, = 
all k G /C. Since for each sine in the sum 

f co \ f 1 w = 2tt/c7t, jf = jfe 

sine - A; = < 

V 27r /^ / [0 w = 2-irie/T, kl ^ k, 

the filter G{oj) satisfies (ITTb by construction. Switching to the time domain 

g{t) = rect (*-) ^ 6^' 2 ^^, 



which is clearly a time compact filter with support r. 
The SoS class in ( fl9l ) may be extended to 

T 



G(u) 



2ir ^ "' \2tt/t 



where bk / 0, k G K, and <j){uj) is any function satisfying: 

1 uo = 



-ft G /C and 6 fc = 6* for 



(18) 



(19) 



(20) 



(21) 



|w|gN 
arbitrary otherwise. 

This more general structure allows for smooth versions of the rect function, which is important when 
practically implementing analog filters. 

The function g(t) represents a class of filters determined by the parameters {frfcjfcg/c. These degrees 
of freedom offer a filter design tool where the free parameters {bk}k<=K ma y be optimized for different 
goals, e.g., parameters which will result in a feasible analog filter. In Theorem [2] below, we show how 
to choose {bk} to minimize the mean-squared error (MSE) in the presence of noise. 
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(a) Periodic signal 

1.5 i 
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(b) Sampling filter 



■ Filter output 
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0.2 0.4 0.6 0.8 

time [units of t] 

(c) Low rate samples 



Fig. 4. Compressed samples of pulse streams (a) Original periodic signal consisting of 5 Gaussians (3 periods are shown), (b) 
Sampling filter, (c) Low rate samples depicted over the filtered signal. 



Determining the parameters {bk}keic m ay be viewed from a more empirical point of view. The impulse 
response of any analog filter having support r may be written in terms of a windowed Fourier series as 

$(t) = rect (-} (3 k e j2nkt/T . (22) 
Confining ourselves to filters which satisfy /3k ^ 0, k G JC, we may truncate the series and choose: 

* * eK (23, 
k i K 

as the parameters of g(t) in ( fT9l ). With this choice, g(t) can be viewed as an approximation to 5>(t). 
Notice that there is an inherent tradeoff here: using more coefficients will result in a better approximation 
of the analog filter, but in turn will require more samples, since the number of samples ./V must be greater 
than the cardinality of the set K. 
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Fig. 5. (a) Estimated time-delays and amplitudes depicted over the original signal, (b) Reconstructed signal vs. original one. 
The reconstruction is exact to numerical precision. 



To demonstrate the filter g(t) we first choose /C = {— p, . . . ,p} and set all coefficients {bk} to one, 
resulting in 

g (t) = rect £ = rect (7) D p (2*t/T), (24) 

where the Dirichlet kernel D p (t) is defined by 

U P^~L, e ~ sin(t/2 ) ■ 

The resulting filter for p = 10 and r = 1 sec, is depicted in Fig. [2 This filter is also optimal in an MSE 
sense for the case h(t) = 6(t), as we show in Theorem |2] In Fig. [3] we plot g(t) for the case in which 
the fefc's are chosen as a length-M symmetric Hamming window: 

b k = 0.54 - 0.46 cos ^ 2?r fc + l^/ 2 J ^ ; k E JC. (26) 

Notice that in both cases the coefficients satisfy bk = b*_ k , and therefore, the resulting filters are real 
valued. 

In the presence of noise, the choice of {bk}keK will effect the performance. Consider the case in which 
digital noise is added to the samples c, so that y = c + w, with w denoting a white Gaussian noise 
vector. Using (fT3l ). 

y = V(-t s )Bx + w (27) 

where B is a diagonal matrix, having {b k } on its diagonal. To choose the optimal B we assume that 
the {ai} are uncorrected with variance a\, independent of {tj}, and that {fy} are uniformly distributed 
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in [0,t). Since the noise is added to the samples after filtering, increasing the filter's amplification will 
always reduce the MSE. Therefore, the filter's energy must be normalized, and we do so by adding the 
constraint Tr(B*B) = 1. Under these assumptions, we have the following theorem: 



Theorem 2. The minimal MSE of a linear estimator of x from the noisy samples y in (1271 ) is achieved 
by choosing the coefficients 

J2 



\h\ 2 



^ A > \h t \ 4 N/a 2 

where = H(2irk/T)o~ a \/~L/T and are arranged in an increasing order of \hf.\, 



(28) 



^ (w-W , (29) 

i=m+l 

and m is the smallest index for which A < |fo m +i| 4 Af/<7 2 . 

Proof: See the Appendix. ■ 
An important consequence of Theorem 2 is the following corollary. 

Corollary 1. If {h^ 2 = \h?\ 2 , Vk,£ £ fC then the optimal coefficients are \bi\ 2 = 1/|/C|, \/k G K.. 



Proof: It is evident from (1281 ) that if \hk\ = \h(\ then \bk\ = \bg\. To satisfy the trace constraint 
Tr(B*B) = 1, A cannot be chosen such that all h = 0. Therefore, \bi\ 2 = l/\K\ for all i G K,. ■ 

From Corollary [T] it follows that when h(t) = 5(t), the optimal choice of coefficients is bk = bj for 
all k and j. We therefore use this choice when simulating noisy settings in the next section. 

Our sampling scheme for the periodic case consists of sampling kernels having compact support in the 
time domain. In the next section we exploit the compact support of our filter, and extend the results to 
the finite stream case. We will show that our sampling and reconstruction scheme offers a numerically 
stable solution, with high noise robustness. 

D. Simulations 

1) Demonstration of Our Sampling Scheme: To demonstrate our results, we consider an input x(t) 
consisting of L = 5 delayed and weighted versions of a Gaussian pulse 

h(t) = -jL= exp(-t 2 /2a 2 ), (30) 
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with parameter a = 7 ■ 10 -3 , and period r = 1. The time-delays and amplitudes were chosen randomly. 
In order to demonstrate near-critical sampling we choose the set of indices K, = {— L,...,L} with 
cardinality M = \K,\ = 11. We filter x(t) with g(t) of d2"6l ). The filter output is sampled uniformly N 
times, with sampling period T = t/N, where N = M = 11. The sampling process is depicted in Fig. |4] 
The vector x is obtained using (fT4l ). and the delays and amplitudes are determined by the annihilating 
filter method. Reconstruction results are depicted in Fig. [5] The estimation and reconstruction are both 
exact to numerical precision. 

Analog filtering operations are carried out by discrete approximations over a fine grid. The analog 
signal and filters are mimicked by high rate digital signals. Since the sampling rate which constructs the 
fine grid is between 2-3 orders of magnitude higher than the final sampling rate T, the simulations reflect 
very well the analog results. 

2) Noisy Case: We now consider the case in which the samples are corrupted by noise. Our signal 
consists of L = 2 pulses h(t) = 5(t). The period was set to r = 1, K = {—2, . . . , 2}, and N = M = 5 
samples were taken, sampled uniformly with sampling period T = t/N. We choose g(t) given by (l24l ). 
As explained earlier, only the values of the filter at points 2irk/T, k G K affect the samples (see (fTTTl). 
Since the values of the filter at the relevant points coincide and are equal to one for the low-pass filter 
Q and g*(—t), the resulting samples for both settings are identical. Therefore, we present results for 
our method only, and state that the exact same results are obtained using the approach of Q. 

In our setup white Gaussian noise (AWGN) with variance a\ is added to the samples, where we define 
the SNR as: 

-Hell 2 

SNR= NU ' , (31) 

with c denoting the clean samples. In our experiments the noise variance is set to give the desired SNR. 

The simulation consists of 1000 experiments for each SNR, where in each experiment a new noise 
vector is created. We choose t = r • (1/3 2/3) T and a = r • (1 l) T , where these vectors remain constant 
throughout the experiments. We define the error in time-delay estimation as the average of ||t — 1 1 1 § , 
where t and t denote the true and estimated time-delays, respectively, sorted in increasing order. The 
error in amplitudes is similarly defined by ||a — a^. In Fig.|6]we show the error as a function of SNR for 
both delay and amplitude estimation. Estimation of the time-delays is the main interest in FRI literature, 
due to special nonlinear methods required for delay recovery. Once the delays are known, the standard 
least-squares method is typically used to recover the amplitudes, therefore, we focus on delay estimation 
in the sequel. 
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SNR [dB] SNR [dB] 

(a) (b) 
Fig. 6. Performance as a function of SNR, using our periodic approach. Estimation error in (a) delays, and (b) amplitudes. 

Finally, for the same setting we can improve reconstruction accuracy at the expense of oversampling, 
as illustrated in Fig. [7] Here we show recovery performance for oversampling factors of 1, 2, 4 and 8. 
The oversampling was exploited using the total least-squares method, followed by Cadzow's iterative 
denoising (both described in detail in ifTOl "). 




10 20 30 40 50 



SNR [dB] 



Fig. 7. The effect of oversampling on estimation error. Oversampling by a factor of 1, 2, 4 and 8. 



16 



III. Finite Stream of Pulses 

A. Extension of SoS Class 

Consider now a finite stream of pulses, denned as 

L 

x(t) = J2 a Mt-ti), tie[0,T),aieR,l = l...L, (32) 
1=1 

where, as in Section [TTJ h(t) is a known pulse shape, and {tua{\f =l are the unknown delays and 
amplitudes. The time-delays {ti}i =1 axe restricted to lie in a finite time interval [0,r). Since there are 
only 2L degrees of freedom, we wish to design a sampling and reconstruction method which perfectly 
reconstructs x(t) from 2L samples. In this section we assume that the pulse h(t) has finite support R, 
i.e., 

h(t) = 0, V|t| > R/2. (33) 

This is a rather weak condition, since our primary interest is in very short pulses which have wide, or 
even infinite, frequency support, and therefore cannot be sampled efficiently using classical sampling 
results for bandlimited signals. We now investigate the structure of the samples taken in the periodic 
case, and design a sampling kernel for the finite setting which obtains precisely the same samples c[n], 
as in the periodic case. 

In the periodic setting, the resulting samples are given by ( fTOb . Using git) of ( fT9l ) as the sampling 
kernel we have 

C N = (9(t - nT),x(t)) 

L poo 

= ^2^2 ai h{t - ti - mr)g* {t - nT)dt 

m£Z 1=1 

L poo 

= ^2Y. ai h(t)g*(t-(nT-ti-mr))dt 

m£Z 1=1 
L 

= ^2^2aiip{nT -t h - mr), (34) 



mSZ 1=1 



where we defined 



<p(#) = {g{t-&),h(t))' (35) 
Since g(t) in (fT9l ) vanishes for all \t\ > r/2 and h(t) satisfies 031 , the support of tp(t) is (R + r), i.e., 

<p{t) = for all \t\ > (R + t)/2. (36) 
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Using this property, the summation in (1341 ) will be over nonzero values for indices m satisfying 

\nT-ti -mr\ < (R + t)/2. (37) 

Sampling within the window [0, r), i.e., nT G [0, r), and noting that the time-delays lie in the interval 
t t G [0,r), I = 1 . . . L, ((37) implies that 

{R + t)/2 > \nT -t h - mr| > |m|r - |nT - t t \ > (|m| - l)r. (38) 

Here we used the triangle inequality and the fact that \nT — ty\ < r in our setting. Therefore, 

'R/t + 3 



, , ^/r + 3 . . 

\m\ < => m < 

ii 2 ii- 



1 = r, (39) 



i.e., the elements of the sum in (1341 ) vanish for all m but the values in (I39I ). Consequently, the infinite 
sum in (l34l reduces to a finite sum over m < |r| so that (l34l) becomes 

r L 

c [ n ] = 5^ o,np{nT — t[ — mr) 

m=—r 1=1 
L 



' /'OO 

XI Xl a W h^-^g*^-^ + mT)dt 



m=—r 1=1 

= /^git-nT + m^^aMt-t^Y (40) 
\m=— r /=1 / 

where in the last equality we used the linearity of the inner product. Defining a function which consists 
of (2r + 1) periods of g(t): 

r 

9r{t) = 9(t + mr), (41) 

m=—r 

we conclude that 

c[n] = (g r (t-nT),x(t)}. (42) 

Therefore, the samples c[n] can be obtained by filtering the aperiodic signal x(t) with the filter g*(—t) 
prior to sampling. This filter has compact support equal to (2r + l)r. Since the finite setting samples 
(|42l are identical to those of the periodic case (l34l . recovery of the delays and amplitudes is performed 
exactly the same as in the periodic setting. 

We summarize this result in the following theorem. 

Theorem 3. Consider the finite stream of pulses given by: 

L 

x(t) = X aih(t - ti), ti e [0, r), ai G R, 
l=i 



18 



where h(t) has finite support R. Choose a set JC of consecutive indices for which H(2irk/r) 7^ 0, VA; G /C. 
Then, N samples given by: 

c [n] = (g r (t -nT),x{t)), n = . . . N - 1, nT e [0, r), 

where r is defined in (l39l ). arccf g r (t) itf compactly supported and defined by (|4TI) (based on the filter g(t) 
in (fTTTl). uniquely determine the signal x(t) as long as N > \fC\ > 2L. 

If, for example, the support of /i(t) satisfies R < t then we obtain from (|39l ) that r = 1. Therefore, 
the filter in this case would consist of 3 periods of g(t): 

93p(t) = g r (t)\ r=1 =g(t-T)+g(t) + g(t + T). (43) 

Practical implementation of the filter may be carried out using delay-lines. The relation of this scheme 
to previous approaches will be investigated in Section [V] 

B. Simulations 

1) Demonstration of the Sampling Scheme: The input signal x(t) consists of L = 5 delayed and 
weighted versions of the pulse h(t) = S(t). The delays and weights were chosen randomly. We choose 
K, = {— L, . . . , L}, so that M = \JC\ = 11. Since the support of h(t) satisfies R < r the parameter r in 
(l39l ) equals 1, and therefore we filter x(t) with gz p {t) defined in (l43l . The coefficients bk, k G K, were 
all set to one. The output of the filter is sampled uniformly N times, with sampling period T = t/N, 
where N = M = 11. Perfect reconstruction is achieved as can be seen in Fig. [8] The estimation is exact 
to numerical precision. 

2) High Order Problems: The same simulation was carried out with L = 20 diracs. The results are 
shown in Fig. [9] Here again, the reconstruction is perfect even for large L. 

3) Noisy Case: We now consider the performance of our method in the presence of noise. In addition, 
we compare our performance to the B -spline and E-spline methods proposed in |[T3l , and to the Gaussian 
sampling kernel Q. We examine 4 scenarios, in which the signal consists of L = 2, 3, 5, 20 dirac|]- In 
our setup, the time-delays are equally distributed in the window [0, r), with r = 1, and remain constant 
throughout the experiments. All amplitudes are set to one. 

'Due to computational complexity of calculating the time-domain expression for high order E-splines, the functions were 
simulated up to order 9, which allows for L = 5 pulses. 
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Fig. 8. Application of the filter ga p (t) on a finite stream of L — 5 diracs. 
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Fig. 9. High order problems: application of the filter ga p (t) on a finite stream of L = 20 diracs. 



The index set of the SoS filter is K, = {— L, . . . , L}. Both B-splines and E-splines are taken of order 
2L — 1, and for E-splines we use purely imaginary exponents, equally distributed around the complex 
unit circle. The sampling period for all methods is T = t/N. 

The method of noise corruption is the same as in Section III-D2I In order to maintain the same SNR 
conditions throughout all methods, the noise level is chosen with respect to the resulting sequence of 
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Fig. 10. Performance in the presence of noise: finite stream case. Our method vs. B-spline, E-spline 1131 and Gaussian t3l 
sampling kernels, (a) L = 2 dirac pulses are present, (b) L = 3 pulses, (c) high value of L — 5 pulses, and (d) the performance 
for a very high value of L = 20 (without E-spline simulation, due to computational complexity of calculating the time-domain 
expression for high values of L). 



samples. In other words, a n in (f3TT > is method-dependent, and is determined by the desired SNR and the 
samples of the specific technique. Hard thresholding was implemented in order to improve the spline 
methods, as suggested by the authors in lfl"3l . The threshold was chosen to be 3a n , where a n is the 
standard deviation of the AWGN. For the Gaussian sampling kernel the parameter a was optimized and 
took on the value of a = 0.25, 0.28, 0.32, 0.9, respectively. 

The results are given in Fig. [10] For L = 2 all methods are stable, where E-splines exhibit better 
performance than B-splines, and Gaussian and SoS approaches demonstrate the lowest errors. As the 
value of L grows, the advantage of the SoS filter becomes more prominent, where for L > 5, the 
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performance of Gaussian and both spline methods deteriorate and have errors approaching the order of 
r. In contrast, the SoS filter retains its performance nearly unchanged even up to L = 20, where the 
B-spline and Gaussian methods are unstable. The improved version of the Gaussian approach presented 
in lfl2l would not perform better in this high order case, since it fails for L > 9, as noted by the authors. 
A comparison of our approach to previous methods will be detailed in Section [V] 

IV. Infinite Stream of Pulses 
We now consider the case of an infinite stream of pulses 

z(t) = J2 a Mt-tl), ti^ieR. (44) 

We assume that the infinite signal has a bursty character, i.e., the signal has two distinct phases: a) bursts 
of maximal duration r containing at most L pulses, and b) quiet phases between bursts. For the sake of 
clarity we begin with the case h(t) = 5{t). For this choice the filter g*(—t) in (|4TT> reduces to gl p (— t) 
of 63. 

Since the filter 5 , L(— t) has compact support 3r we are assured that the current burst cannot influence 
samples taken 3r/2 seconds before or after it. In the finite case we have confined ourselves to sampling 
within the interval [0, r). Similarly, here, we assume that the samples are taken during the burst duration. 
Therefore, if the minimal spacing between any two consecutive bursts is 3r/2, then we are guaranteed that 
each sample taken during the burst is influenced by one burst only, as depicted in Fig. [TT] Consequently, 
the infinite problem can be reduced to a sequential solution of local distinct finite order problems, 
as in Section |III] Here the compact support of our filter comes into play, allowing us to apply local 
reconstruction methods. 









93 P {t) filter support = 3r 






























i 


. 1 


t 


I 


t 


; ) 

i -T "4 


k 


It 


■• • • 

— ^ — 



-0.hr t Ut 2.5r 3.5r 



1st burst 2nd burst 

Fig. 11. Bursty signal z(t). Spacing of St/2 between bursts ensures that the influence of the current burst ends before taking 
the samples of the next burst. This is due to the finite support, 3r of the sampling kernel g*j P {—t). 
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In the above argument we assume we know the locations of the bursts, since we must acquire samples 
from within the burst duration. Samples outside the burst duration are contaminated by energy from 
adjacent bursts. Nonetheless, knowledge of burst locations is available in many applications such as 
synchronized communication where the receiver knows when to expect the bursts, or in radar or imaging 
scenarios where the transmitter is itself the receiver. 

We now state this result in a theorem. 

Theorem 4. Consider a signal z(t) which is a stream of bursts consisting of delayed and weighted 
diracs. The maximal burst duration is r, and the maximal number of pulses within each burst is L. Then, 
the samples given by 

C M = (93p(t ~ nT ), z (t))i n£Z 

where g3 P (t) is defined by (143b . are a sufficient characterization of z(t) as long as the spacing between 
two adjacent bursts is greater than 3r/2, and the burst locations are known. 

Extending this result to a general pulse h(t) is quite straightforward, as long as h(t) is compactly 
supported with support R, and we filter with g*(—t) as defined in (l4lT ) with the appropriate r from j39T l. 
If we can choose a set /C of consecutive indices for which H(27rk/r) / 0, Vfc 6 K, and we are guaranteed 
that the minimal spacing between two adjacent bursts is greater than ((2r + l)r + R) /2, then the above 
theorem holds. 

V. Related Work 

In this section we explore the relationship between our approach and previously developed solutions 

a, ma, ma, m. 

A. Periodic Case 

The work in (3l was the first to address efficient sampling of pulse streams, e.g., diracs. Their approach 
for solving the periodic case was ideal lowpass filtering, followed by uniform sampling, which allowed to 
obtain the Fourier series coefficients of the signal. These coefficients are then processed by the annihilating 
filter to obtain the unknown time-delays and amplitudes. In Section [TTJ we derived a general condition on 
the sampling kernel (fTTb . under which recovery is guaranteed. The lowpass filter of (3] is a special case 
of this result. The noise robustness of both the lowpass approach and our more general method is high 
as long as the pulses are well separated, since reconstruction from Fourier series coefficients is stable in 
this case. Both approaches achieve the minimal number of samples. 
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The lowpass filter is bandlimited and consequently has infinite time-support. Therefore, this sampling 
scheme is unsuitable for finite and infinite streams of pulses. The SoS class introduced in Section [HI 
consists of compactly supported filters which is crucial to enable the extension of our results to finite 
and infinite streams of pulses. A comparison between the two methods is shown in Table U 

TABLE I 

Periodic case - Comparison with previous work 



Feature 


Lowpass filter 
® 


Proposed method 


Degrees of freedom 


2L 


No. of samples 


2L + 1 


2L 


Time- support 


Infinite 


t, finite support 
allows extension 
to finite & infinite 
cases 


Noise Robustness 


High 


High 


Analog implementa- 
tion 


Approximate 
lowpass filter 


Approximate 
finite support 
filter 



B. Finite Pulse Stream 

The authors of O proposed a Gaussian sampling kernel for sampling finite streams of Diracs. The 
Gaussian method is numerically unstable, as mentioned in lfl2l . since the samples are multiplied by a 
rapidly diverging or decaying exponent. Therefore, this approach is unsuitable for L > 6. Modifications 
proposed in |[T2l exhibit better performance and stability. However, these methods require substantial 
oversampling, and still exhibit instability for L > 9. 

In |fT3l the family of polynomial reproducing kernels was introduced as sampling filters for the model 
d32l ). B-splines were proposed as a specific example. The B-spline sampling filter enables obtaining 
moments of the signal, rather than Fourier coefficients. The moments are then processed with the same 
annihilating filter used in previous methods. However, as mentioned by the authors, this approach is 
unstable for high values of L. This is due to the fact that in contrast to the estimation of Fourier 
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coefficients, estimating high order moments is unstable, since unstable weighting of the samples is earned 
out during the process. 

Another general family introduced in |fT3l for the finite model is the class of exponential reproducing 
kernels. As a specific case, the authors propose E-spline sampling kernels. The CTFT of an E-spline of 
order N + 1 is described by 

N 

P a (u) = n , (45) 

where a = (a>o, a± , . . . , ayv) are free parameters. In order to use E-splines as sampling kernels for pulse 
streams, the authors propose a specific structure on the a's, a n = «o + Choosing exponents having 
a non-vanishing real part results in unstable weighting, as in the B-spline case. However, choosing the 
special case of pure imaginary exponents in the E-splines, already suggested by the authors, results in 
a reconstruction method based on Fourier coefficients, which demonstrates an interesting relation to our 
method. The Fourier coefficients are obtained by applying a matrix consisting of the exponent spanning 
coefficients {c m)n }, (see |fT3l ), instead of our Vandermonde matrix relation (fT4l ). With this specific choice 
of parameters the E-spline function satisfies (fTTT) . 

Interestingly, with a proper choice of spanning coefficients, it can be shown that the SoS class can 
reproduce exponentials with frequencies \2irk / t}i- & jc, and therefore satisfies the general exponential 
reproduction property of |[T3l . However, the SoS filter proposes a new sampling scheme which has 
substantial advantages over existing methods including E-splines. The first advantage is in the presence 
of noise, where both methods have the following structure: 

y = Ax + w, (46) 

where w is the noise vector. While the Fourier coefficients vector x is common to both approaches, 
the linear transformation A is method dependent, and therefore the sample vector y is different. In our 
approach with g(t) of (1241 . A is the DFT matrix, which for any order L has a condition number of 
1. However, in the case of E-splines the transformation matrix A consists of the E-spline exponential 
spanning coefficients, which has a much higher condition number, e.g., above 100 for L = 5. Conse- 
quently, some Fourier coefficients will have much higher values of noise than others. This scenario of 
high variance between noise levels of the samples is known to deteriorate the performance of spectral 
analysis methods ifTTll . the annihilating filter being one of them. This explains our simulations which 
show that the SoS filter outperforms the E-spline approach in the presence of noise. 

When the E-spline coefficients a axe pure imaginary, it can be easily shown that ( [431 ) becomes a 
multiplication of shifted sines. This is in contrast to the SoS filter which consists of a sum of sines in 
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the frequency domain. Since multiplication in the frequency domain translates to convolution in the time 
domain, it is clear that the support of the E-spline grows with its order, and in turn with the order of 
the problem L. In contrast, the support of the SoS filter remains unchanged. This observation becomes 
important when examining the infinite case. The constraint on the signal in |[T3l is that no more than L 
pulses be in any interval of length LPT, P being the support of the filter, and T the sampling period. Since 
P grows linearly with L, the constraint cast on the infinite stream becomes more stringent, quadratically 
with L. On the other hand, the constraint on the infinite stream using the SoS filter is independent of L. 

We showed in simulations that typically for L > 5 the estimation errors, using both B-spline and E- 
spline sampling kernels, become very large. In contrast, our approach leads to stable reconstruction even 
for very high values of L, e.g., L = 100. In addition, even for low values of L we showed in simulations 
that although the E-spline method has improved performance over B-splines, the SoS reconstruction 
method outperforms both spline approaches. A comparison is described in Table HO 

TABLE II 
Finite case - comparison 



Feature 


Gaussian 
filter (3) 


Spline Fil- 
ter fT3l 


Proposed 
method 


Degrees of free- 
dom 


2L 


No. of samples 


2L 


Time-support 


Infinite 


Finite 


Finite 


Stability 


Unstable 
for L > 6 


Unstable 
for L > 5 


Stable even 
for L = 
100 


Noise Robustness 


Low 


Low 


High 



C. Infinite Streams 

The work in |[T3l addressed the infinite stream case, with h(t) = 5(t). They proposed filtering the 
signal with a polynomial reproducing sampling kernel prior to sampling. If the signal has at most L 
diracs within any interval of duration LPT, where P denotes the support of the sampling filter and T the 
sampling period, then the samples are a sufficient characterization of the signal. This condition allows 
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to divide the infinite stream into a sequence of finite case problems. In our approach the quiet phases of 
1.5r between the bursts of length r enable the reduction to the finite case. 

Since the infinite solution is based on the finite one, our method is advantageous in terms of stability 
in high order problems and noise robustness. However, we do have an additional requirement of quiet 
phases between the bursts. 

Regarding the sampling rate, the number of degrees of freedom of the signal per unit time, also known 
as the rate of innovation, is p = 2L/2.5r, which is the critical sampling rate. Our sampling rate is 2L/t 
and therefore we oversample by a factor of 2.5. In the same scenario, the method in |fT3l would require 
a sampling rate of LP '/2.5r, i.e., oversampling by a factor of P/2. Properties of polynomial reproducing 
kernels imply that P > 2L, therefore for any L > 3, our method exhibits more efficient sampling. A 
table comparing the various features is shown in Table |III] 

Recent work lfl4l presented a low complexity method for reconstructing streams of pulses (both infinite 
and finite cases) consisting of diracs. However the basic assumption of this method is that there is at 
most one dirac per sampling period. This means we must have prior knowledge about a lower limit on 
the spacing between two consecutive deltas, in order to guarantee correct reconstruction. In some cases 
such a limit may not exist; even if it does it will usually force us to sample at a much higher rate than 
the critical one. 



TABLE III 
Infinite case - Comparison 



Feature 


Spline filter El 


Proposed method 


Signal model 


No more than L 
pulses in any inter- 
val of LPT sec 


Bursty character: 
burst - r, quiet 
phase 1.5r 


Rate of innovation 


p = 2L/2.5r 


Sampling rate 


P- p/2 


2.5p 


For L > 3 =>■ 
P/2 > 3 


Proposed sampling scheme is more efficient 


Noise Robustness 


Low 


High 


Stability 


Unstable for L > 

5 


Stable for L = 
100 
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VI. Application - Ultrasound Imaging 

An interesting application of our framework is ultrasound imaging. In ultrasonic imaging an acoustic 
pulse is transmitted into the scanned tissue. The pulse is reflected due to changes in acoustic impedance 
which occur, for example, at the boundaries between two different tissues. At the receiver, the echoes are 
recorded, where the time-of-arrival and power of the echo indicate the scatterer's location and strength, 
respectively. Accurate estimation of tissue boundaries and scatterer locations allows for reliable detection 
of certain illnesses, and is therefore of major clinical importance. The location of the boundaries is often 
more important than the power of the reflection. This stream of pulses is finite since the pulse energy 
decays within the tissue. We now demonstrate our method on real 1-dimensional (ID) ultrasound data. 

The multiple echo signal which is recorded at the receiver can be modeled as a finite stream of pulses, 
as in d32l) . The unknown time-delays correspond to the locations of the various scatterers, whereas the 
amplitudes correspond to their reflection coefficients. The pulse shape in this case is a Gaussian defined 
in OUT ), due the physical characteristics of the electro-acoustic transducer (mechanical damping). We 
assume the received pulse-shape is known, either by assuming it is unchanged through propagation, 
through physically modeling ultrasonic wave propagation, or by prior estimation of received pulse. Full 
investigation of mismatch in the pulse shape is left for future research. 

In our setting, a phantom consisting of uniformly spaced pins, mimicking point scatterers, was scanned 
by GE Healthcare's Vivid-i portable ultrasound imaging system |[20l . l2Tl . using a 3S-RS probe. We use 
the data recorded by a single element in the probe, which is modeled as a ID stream of pulses. The 
center frequency of the probe is f c = 1.7021 MHz, The width of the transmitted Gaussian pulse in this 
case is a = 3 • 10~ 7 sec, and the depth of imaging is R max = 0.16 m corresponding to a time window 
ofl r = 2.08 • 10~ 4 sec. 

In this experiment all filtering and sampling operations are carried out digitally in simulation. The 
analog filter required by the sampling scheme is replaced by a lengthy Finite Impulse Response (FIR) 
filter. Since the sampling frequency of the element in the system is f s = 20 MHz, which is more 
than 5 times higher than the Nyquist rate, the recorded data represents the continuous signal reliably. 
Consequently, digital filtering of the high-rate sampled data vector (4160 samples) followed by proper 
decimation mimics the original analog sampling scheme with high accuracy. The recorded signal is 
depicted in Fig. |T2j The band-pass ultrasonic signal is demodulated to base-band, i.e., envelope-detection 
is performed, before inserted into the process. 



2 The speed of sound within the tissue is 1550 m/sec. 
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Fig. 12. Recorded ultrasound imaging signal. The data was acquired by GE healthcare's Vivid-i ultrasound imaging system. 

We carried out our sampling and reconstruction scheme on the aforementioned data. We set L = 4, 
looking for the strongest 4 echoes. Since the data is corrupted by strong noise we over-sampled the 
signal, obtaining twice the minimal number of samples. In addition, hard-thresholding of the samples 
was implemented, where we set the threshold to 10 percent of the maximal value. We obtained N = 17 
samples by decimating the output of the lengthy FIR digital filter imitating <?L(— t) from d43l ), where 
the coefficients {b^} were all set to one. In Fig. [T3h the reconstructed signal is depicted vs. the full 
demodulated signal using all 4160 samples. Clearly, the time-delays were estimated with high precision. 
The amplitudes were estimated as well, however the amplitude of the second pulse has a large error. 
This is probably due to the large values of noise present in its vicinity. However, as mentioned earlier, 
the exact locations of the scatterers is often more important than the accurate reflection coefficients. We 
carried out the same experiment only now oversampling by a factor of 4, resulting in N = 33 samples. 
Here no hard-thresholding is required. The results are depicted in Fig. [T"3b . and are very similar to our 
previous results. In both simulations, the estimation error in the pulse location is around 0.1 mm. 

Current ultrasound imaging technology operates at the high rate sampled data, e.g., f s = 20 MHz in 
our setting. Since there are usually 100 different elements in a single ultrasonic probe each sampled at 
a very high rate, data throughput becomes very high, and imposes high computational complexity to the 
system, limiting its capabilities. Therefore, there is a demand for lowering the sampling rate, which in 
turn will reduce the complexity of reconstruction. Exploiting the parametric point of view, our sampling 
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Fig. 13. Applying our gs P (t) filter method on real ultrasound imaging data. Results are shown vs. full demodulated signal 
which uses all 4160 samples. Reconstructed signal (a) using N — 17 samples only and hard-thresholding, and (b) using N — 33 
samples without thresholding. 



scheme reduces the sampling rate by 2 orders of magnitude, from 4160 to around 30 samples in our 
setting, while estimating the locations of the scatterers with high accuracy. 

VII. Conclusions 

We presented efficient sampling and reconstruction schemes for streams of pulses. For the case of 
a periodic stream of pulses, we derived a general condition on the sampling kernel which allows a 
single-channel uniform sampling scheme. Previous work (3] is a special case of this general result. We 
then proposed a class of filters, satisfying the condition, with compact support. Exploiting the compact 
support of the filters, we constructed a new sampling scheme for the case of a finite stream of pulses. 
Simulations show this method exhibits better performance than previous techniques 0, |fT3l , in terms 
of stability in high order problems, and noise robustness. An extension to an infinite stream of pulses 
was also presented. The compact support of the filter allows for local reconstruction, and thus lowers 
the complexity of the problem. Finally, we demonstrated the advantage of our approach in reducing the 
sampling and processing rate of ultrasound imaging, by applying our techniques to real ultrasound data. 
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Appendix 
Proof of Theorem [2] 

The MSE of the optimal linear estimator of the vector x from the measurement vector y is known to 
be El 

MSE = Tr {K ,.,.} - Tr { R J R ;/ ,. } . (47) 
The covariance matrices in our case are 

R-xy = R-ia;B*V* (48) 

R yy = VBR ra B* V* + a 2 I, (49) 

where we used (l27l ). and the fact that H ww = a 2 \ since w is a white Gaussian noise vector. Under our 
assumptions on {ti} and {ai}, denoting h k = Hfiirk/r), and using ((5]) 



(R xx ) k>k , = E{X[k]X*[k'}} 



L L 



1=1 i'=i 

lh k h k ,Y J E{e-^ k - k '^} 



T 

1=1 



[ T -e-^ k - k ^dt 

T Jo T 



i=i 

2 



°-±L\h k \ 2 b Kk ,. (50) 



Denoting by H a diagonal matrix with Mi element |/ifc| 2 = \h k \ 2 a 2 L/T 2 we have 

R-xx = H. (51) 

Since the first term of (1471 is independent of B, minimizing the MSE with respect to B is equivalent 
to maximizing the second term in d47l ). Substituting (|48l,d49l and ( [511 ) into this term, the optimal B is 
a solution to 

maxTr{HB*V*(VBHB*V* +cr 2 I)^ 1 VBH} (52) 
s.t. Tr(B*B) = 1. 
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Using the matrix inversion formula E3l . 



(VBHB*V* + a 2 !)' 1 
= — ( I — VB ( a 2 ^ 1 + B*V*Vb) 



It is easy to verify from the definition of V in (fT3l that 



N-l 



(53) 



(54) 



1=0 



Therefore, the objective in ( 1521 ) equals 



<! ^ HB* | I - B ( ^-H 1 + B*B ) 



£n 2 (i 

i=l v 



N 
a 2 /N 



(55) 



\bi\ 2 \hi\ 2 + a 2 /N , 

where we used the fact that B and H are diagonal. 

We can now find the optimal B by maximizing (|55T ). which is equivalent to minimizing the negative 
term: 



\K\ 



\hA 2 

mm 

B \h-\2\h-\2M/rr2 



\K\ 

S.t. Y^\ b i 



Denoting fa = \bi\ 2 , 



subject to 



^l + \bi\ 2 \hi\ 2 N/a 
becomes a convex optimization problem: 

^l + fafh^N/a 2 



fa >o 



1. 



(56) 



mm 



|/C| 
i=l 

To solve (1571 ) subject to (158T ) and d59l ), we form the Lagrangian: 



(57) 

(58) 
(59) 



\K\ \ \K\ 



(60) 



where from the Karush-Kuhn-Tucker (KKT) conditions ll24l . > and = 0. Differentiating ([60 
with respect to and equating to 

\hi\ A N/a 2 



(l + fa\hi\ 2 N/a 2 ) 2 



+ Vi = A, 



(61) 
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so that A > 0, since hi > by construction of H (see Theorem [T]). If A > \hi\ A N/a 2 then in > 0, and 
therefore, ft = from KKT. If A < \hi\ 4 N/a 2 then from <gT) = and 



■2 




1 



12 



(62) 



The optimal ft is therefore 



iV 



1 ' A < N 4 iV/a 2 

(63) 



IM 2 , 

A > \hi\ 4 N/a 2 

where A > is chosen to satisfy d59l ). Note that from (l63l . if ft 7^ and i < j, then ft- / as well, 
since |/tj| axe in an increasing order. We now show that there is a unique A that satisfies d59l ). Define the 
function 

\K\ 

S(A) = £>(A) - 1, (64) 

8=1 

so that A is a root of Q{\). Since the |/tj|'s are in an increasing order, |/j|;q| = maxj It is clear 
from (l63l that ^(A) is monotonically decreasing for < A < | 4 A^/cr 2 . In addition, <5(A) = —1 for 
A > |/t|,q| 4 iV/o- 2 , and ^(A) > for A ->■ 0. Thus, there is a unique A for which ( f59b is satisfied. 

Substituting (l63l) into d59l) , and denoting by m the smallest index for which A < |/i m+ i| 4 '/ a 2 , we 
have 



✓A= , (65) 



2 



iV/a 2 + ^ Vl^i 



i=m+l 



completing the proof of the theorem. 
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